This analysis locally maps breakpoint density across the genome by binning breakpoint SV and CNV data. This notebook returns data that are passed on to 02-a-plot-chr-instability-heatmaps.Rmd and 02b-plot-chr-instability-by-histology.Rmd for plotting.
The breaks_density function bins the genome and counts how many chromosome breaks occur for each bin, given the given bin_size.
We will calculate breakpoint density across bins of the genome for each sample and histology group for all three breakpoint datasets set up by 00-setup-breakpoint_data.R.
intersection_of_breaks contains the intersection break counts for both SV and CNV break data.cnv_breaks contains the number of break counts for CNV.sv_breaks contains the number of break counts for SV.This notebook can be run via the command line from the top directory of the repository as follows:
Rscript -e "rmarkdown::render('analyses/chromosomal-instability/01-localization-of-breakpoints.Rmd',
clean = TRUE)"# Change here and it will change the rest
bin_size <- 1e6
# Set seed so heatmaps turn out the same
set.seed(2020)
# Magrittr pipe
`%>%` <- dplyr::`%>%`Read in the custom functions needed for this analysis.
source(file.path("util", "chr-break-calculate.R"))
source(file.path("util", "chr-break-plot.R"))# Path to data directory
data_dir <- file.path("..", "..", "data")
scratch_dir <- file.path("..", "..", "scratch")
# The output from this notebook is used for downstream notebooks so we will save
# it in a special directory
output_dir <- "breakpoint-data"
# Create the output_dir if it does not exist
if (!dir.exists(output_dir)) {
dir.create(output_dir)
}Here’s all the input files we will need:
chr_file <- file.path(data_dir, "WGS.hg38.strelka2.unpadded.bed")
metadata_file <- file.path(data_dir, "pbta-histologies.tsv")
breaks_file <- file.path(output_dir, "breaks_lists.RDS")
uncallable_file <- file.path(
"..", "copy_number_consensus_call", "ref",
"cnv_excluded_regions.bed"
)The histology output file:
output_file <- file.path(output_dir, "histology_breakpoint_binned_counts.RDS")Wrapper function to run break_density for all three datasets: CNV, SV, and intersection data.
all_break_density <- function(breaks_list, sample_ids, return_vector) {
# For given sample(s) ids, run breaks_density for all dataframes items in breaks_list.
#
# Args:
# breaks_list: a data.frame with chromosomal breaks.
# sample_ids: a single or multiple samples that correspond to the sample_col
# argument and will be used for calculating density.
# return_vector: If we don't need the full GenomicRanges object, only the total_counts, put TRUE
#
# Returns:
# GenomicRanges object that contains the densities mapped across the genome
# for the given samples and for all data.frames in breaks_list.
lapply(breaks_list, break_density, # Either intersection, CNV, SV data depending on the iteration.
sample_id = sample_ids,
start_col = "coord",
end_col = "coord",
window_size = bin_size, # Bin size to calculate breaks density
chr_sizes_vector = chr_sizes_vector, # This is the sizes of chromosomes used for binning
unsurveyed_bed = uncallable_bed, # This is the BED file that notes what regions are uncallable
perc_cutoff = .75, # What percentage of each bin needs to be callable for it not to be NA
return_vector = return_vector # If we don't need the full GenomicRanges object
)
}Import data that was set up in breakpoint format in 00-setup-breakpoint-data.R.
breaks_list <- readr::read_rds(breaks_file)Get the list of biospecimen IDs that are in this set.
common_samples <- unique(breaks_list$intersection_of_breaks$samples)Set up metadata
# Read in the metadata
metadata <- readr::read_tsv(metadata_file) %>%
# Only keep metadata for the samples we are working with here
dplyr::filter(Kids_First_Biospecimen_ID %in% common_samples)Parsed with column specification:
cols(
.default = col_character(),
OS_days = [32mcol_double()[39m,
age_last_update_days = [32mcol_double()[39m,
normal_fraction = [32mcol_double()[39m,
tumor_fraction = [32mcol_double()[39m,
tumor_ploidy = [32mcol_double()[39m,
molecular_subtype = [33mcol_logical()[39m
)
See spec(...) for full column specifications.
493 parsing failures.
row col expected actual file
2334 molecular_subtype 1/0/T/F/TRUE/FALSE Group3 '../../data/pbta-histologies.tsv'
2335 molecular_subtype 1/0/T/F/TRUE/FALSE Group4 '../../data/pbta-histologies.tsv'
2336 molecular_subtype 1/0/T/F/TRUE/FALSE Group3 '../../data/pbta-histologies.tsv'
2337 molecular_subtype 1/0/T/F/TRUE/FALSE Group3 '../../data/pbta-histologies.tsv'
2338 molecular_subtype 1/0/T/F/TRUE/FALSE Group3 '../../data/pbta-histologies.tsv'
.... ................. .................. ...... .................................
See problems(...) for more details.Set up chromosome size data that will be used for creating genome bins. It just so happens that this BED file: WGS.hg38.strelka2.unpadded.bed is actually just a list of the chromosome sizes so we are using that for now.
# Set up Chr sizes
chr_sizes_bed <- readr::read_tsv(chr_file,
col_names = c("chrom", "start", "end")
) %>%
# Reformat the chromosome variable to drop the "chr"
dplyr::mutate(chrom = factor(gsub("chr", "", chrom),
levels = c(1:22, "X", "Y", "M")
)) %>%
# Remove sex chromosomes
dplyr::filter(!(chrom %in% c("X", "Y", "M")))Parsed with column specification:
cols(
chrom = [31mcol_character()[39m,
start = [32mcol_double()[39m,
end = [32mcol_double()[39m
)# Make chromosome size named vector for Heatmap annotation
chr_sizes_vector <- chr_sizes_bed$end
names(chr_sizes_vector) <- chr_sizes_bed$chromSet up uncallable regions BED. We will use this to declare NAs in regions with too high of a percentage of uncallable regions.
uncallable_bed <- readr::read_tsv(uncallable_file,
col_names = c("chrom", "start", "end")
) %>%
# Reformat the chromosome variable to drop the "chr"
dplyr::mutate(chrom = factor(gsub("chr", "", chrom),
levels = c(1:22, "X", "Y")
)) %>%
dplyr::filter(
!is.na(chrom),
# Drop sex chromosomes
!(chrom %in% c("X", "Y", "M"))
)Parsed with column specification:
cols(
chrom = [31mcol_character()[39m,
start = [32mcol_double()[39m,
end = [32mcol_double()[39m
)Let’s reformat that the intersection_or_breaks data.frame. It contains both SV data and CNV data for its intersections. We’ll go by SV data here.
# For simplicity's sake, change the name of the intersection column we want to use
breaks_list$intersection_of_breaks <- breaks_list$intersection_of_breaks %>%
dplyr::rename(
coord = sv_ranges.start
)Calculate sample breakpoint densities across the binned genome.
# Get a big list of break densities for each sample.
sample_densities <- lapply(common_samples, function(sample_id) {
all_break_density(
breaks_list = breaks_list,
sample_id = sample_id,
return_vector = TRUE
) # We only want the total counts,
# not the whole GenomicRanges object
})
# Bring along the sample IDs
names(sample_densities) <- common_samplesSwitch list to be by breaks first, samples second.
sample_densities <- purrr::transpose(sample_densities)
# Extract the chromosome names for each bin
chr_bin_names <- names(sample_densities[[1]][[1]])Turn each dataset into its own data.frame and write it to its own TSV file. This will be used by 02-a-plot-chr-instability-heatmaps.Rmd for heatmaps.
# Write the break densities each as their own files
purrr::imap(sample_densities, function(.x, name = .y) {
dplyr::bind_rows(.x) %>%
# Add on the chromosome names so they don't get lost
dplyr::mutate(chr_bin_names) %>%
# Write each to a TSV file
readr::write_tsv(file.path(
output_dir,
paste0(name, "_binned_counts.tsv")
))
})$intersection_of_breaks
$cnv_breaks